#pragma once
#include <vector>
#include <iostream>
#include <lapacke.h>
#include "function.h"
#include "MatrixV.h"
using namespace std;
class twostage_collocation_method
{
private:
    /// 几个参数的关系为h = 1/(m+1), r = (k*nu)/(h*h)
    int m;
    double r, h, nu, k;
    vector<vector<double>> A, I;
public:
    twostage_collocation_method(double _r, double _h) : r(_r), h(_h)
    {
	m = (int)1/h - 1;
	nu = 1;
	k = h*h*r/nu;
	MatrixV M(r, nu, h);
	A = M.A();
	I = M.I();
    };
    /// 根据初始条件函数创建 U^0
    vector<double> U0()
    {
	vector<double> U0;
	for (int i = 0; i < m; i++)
	{
	    U0.push_back(u0((i+1)*h));
	}
	return U0;
    };
    /// 求解器,输入时间步数,保存求解结果到 U
    vector<double> Solver(int t)
    {
	vector<double> MU(2*m), UN(m);
	vector<vector<double>> M1(2*m), U;
	U.push_back(U0());
	vector<vector<double>> A1,A2,A3,A4;
	A1 = -3/4.0*k*A;
	A2 = I + A1;
	A4 = k/12.0*A;
	A3 = I - 5*A4;
	    
	for (int i = 0; i < m; i++)
	{
	    M1[i].resize(2*m);
	    M1[i+m].resize(2*m);
	    for (int j = 0; j < m; j++)
	    {
		M1[i][j] = A1[i][j];
		M1[i][j+m] = A2[i][j];
		M1[i+m][j] = A3[i][j];
		M1[i+m][i+m] = A4[i][j];
	    }
	}
	for (int l = 0; l < t; l++)
	{
	    for (int i = 0; i < m; i++)
	    {
		MU[i] = U[l][i];
		MU[i+m] = U[l][i];
	    }
	    MU = solveU(M1, MU);
	    for (int i = 0; i < m; i ++)
	    {
		UN[i] = MU[i+m];
	    }
	    U.push_back(UN);
	}
	return U[t];
    };
    vector<double> Solver0(int t)
    {
	vector<double> UN(m), MU;
	vector<vector<double>> U, dU, M1, M2;
	U.push_back(U0());
	//dUN = A*U[0];
	dU.push_back(A*U[0]);
	M1 = I - k/2*A;
	M2 = I + k/2*A;
	MU = M2*U[0];
	// double *A1 = new double[m*m];
	// double *F = new double[m];
	// for (int j = 0; j < m; j++)
	// {
	//     F[j] = MU[j];
	//     for (int l = 0; l < m; l++)
	//     {
	// 	A1[j+l*m] = M1[j][l];
	//     }
	// }
	// lapack_int info,m0,n0,lda,ldb,nrhs;
	// m0 = m;
	// n0 = m;
	// nrhs = 1;
	// lda = m;
	// ldb = m;
	// info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',m0,n0,nrhs,A1,lda,F,ldb);
	// for(int i = 0; i < m; i++)
	// {
	//     UN[i] = F[i];
	// }
	UN = solveU(M1,MU);
	U.push_back(UN);
	dU.push_back(1.0/3*A*(4*U[1] - U[0]));
	for (int i = 1; i < t; i++)
	{
	    M1 = I - k/4*A;
	    MU = U[i] + 3/4*k*dU[i];
	    // for (int j = 0; j < m; j++)
	    // {
	    // 	F[j] = MU[j];
	    // 	for (int l = 0; l < m; l++)
	    // 	{
	    // 	    A1[j+l*m] = M1[j][l];
	    // 	}
	    // }
	    // lapack_int info,m0,n0,lda,ldb,nrhs;
	    // m0 = m;
	    // n0 = m;
	    // nrhs = 1;
	    // lda = m;
	    // ldb = m;
	    // info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',m0,n0,nrhs,A1,lda,F,ldb);
	    // for(int i = 0; i < m; i++)
	    // {
	    // 	UN[i] = F[i];
	    // }
	    UN = solveU(M1, MU);
	    U.push_back(UN);
	    dU.push_back(3.0/2*A*U[i+1] - 1.0/2*dU[i]);
	}
	return U[t];
    };
};

class onestage_GaussLegendre_RK_method
{
private:
    /// 几个参数的关系为h = 1/(m+1), r = (k*nu)/(h*h)
    int m;
    double r, h, nu, k;
    vector<vector<double>> A, I;
public:
    onestage_GaussLegendre_RK_method(double _r, double _h) : r(_r), h(_h)
    {
	m = (int)1/h - 1;
	nu = 1;
	k = h*h*r/nu;
	MatrixV M(r, nu, h);
	A = M.A();
	I = M.I();
    };
    /// 根据初始条件函数创建 U^0
    vector<double> U0()
    {
	vector<double> U0;
	for (int i = 0; i < m; i++)
	{
	    U0.push_back(u0((i+1)*h));
	}
	return U0;
    };
    /// 求解器,输入时间步数,保存求解结果到 U
    vector<double> Solver(int t)
    {
	vector<double> MU, UN(m);
	vector<vector<double>> U, dU, M1;
	U.push_back(U0());
	for (int i = 0; i < t; i++)
	{
	    //print(U[i]);
	    //U.push_back((k*A+I)*U[i]);
	    //dU.push_back(dU[i]);
	    MU = U[i];
	    M1 = I - k/2*A;
	    // double *A1 = new double[m*m];
	    // double *F = new double[m];
	    // //cout << "A1,F" << endl;
	    // for (int j = 0; j < m; j++)
	    // {
	    // 	F[j] = MU[j];
	    // 	for (int l = 0; l < m; l++)
	    // 	{
	    // 	    A1[j+l*m] = M1[j][l];
	    // 	}
	    // }
	    // //cout << "A1,F" << endl;
	    // lapack_int info,m0,n0,lda,ldb,nrhs;
	    // m0 = m;
	    // n0 = m;
	    // nrhs = 1;
	    // lda = m;
	    // ldb = m;
	    // info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',m0,n0,nrhs,A1,lda,F,ldb);
	    // //cout << "LAPACK" << endl;
	    // for(int p = 0; p < m; p++)
	    // {
	    // 	UN[p] = F[p];
	    // }
	    //cout << "UN" << endl;
	    UN = solveU(M1,MU);
	    dU.push_back(UN);
	    //print(dU[i]);
	    U.push_back(k*A*dU[i] + U[i]);
	}
	return U[t];
    };
};
